A diffusion Monte Carlo study of small para-Hydrogen clusters 
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Ground state energies and chemical potentials of parahydrogen clusters are calculated from 3 to 40 
molecules using the diffusion Monte Carlo technique with two different P-H2-P-H2 interactions. This 
calculation improves a previous one by the inclusion of three-body correlations in the importance 
sampling, by the time step adjustement and by a better estimation of the statistical errors. Apart 
from the cluster with f3 molecules, no other magic clusters are predicted, in contrast with path 
integral Monte Carlo results. 
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Theoretical studies of parahydrogen clusters have at- 
tracted a growing interest in the past years, partly mo- 
tivated by a recent experiment [l| in which Raman scat- 
tering was used in cryogenic free jets of the pure gas. 
Small changes in the frequency near the Qi{0) line of the 
monomer, were observed and interpreted as intermolecu- 
lar effects on the intramolecular potential. (p-H2)Ar clus- 
ters with = 2 — 8 were clearly identified through fre- 
quency shifts ranging from Av = — 0.40cm~^ for N = 2 
to Ai^ — — 2.35cm~-'^ for iV = 8. The experiment also 
showed a bump at iV = 13, = 33 and N = 55, which 
were interpreted as a signal of magical clusters. However, 
it is worth mentioning that these three values are actu- 
ally extrapolations from smaller clusters and presumably 
are approximate. 

Magical numbers appear in classical Lennard-Jones 
clusters, related to geometrical shapes [2|. Several papers 
have appeared in the last year with the main objective of 
checking the magical numbers found in Ref. 1], and/or 
studying possible superfluidity effects in parahydrogen 
clusters. Indeed, Path Integral Monte Carlo (PIMC) cal- 
culations 3 have found a large superfluid fraction in 
clusters with N=13 and 18 at temperatures T < 2 K. 
A superfluid response has been observed in small clus- 
ters consisting of a carbonyl sulfide cromophore sur- 
rounded by 15-17 P-H2 molecules, all within a large he- 
lium droplet i4J. This has been confirmed by several MC 
simulations [all> 0, B 9 of doped P-H2 clusters. 

Systematic studies of (p-H2)Ar clusters, covering the 
range from = 3 to = 50 molecules, have been 
done based on powerful many-body techniques, as di- 
fusion Monte Carlo (DMC) [10|, PIMC and 
PIMC adapted to the ground state (PIGS)~[lir Whereas 
up to A^ ~ 22 all these calculations are substantially in 
agreement, for heavier clusters there are noticeable dif- 
ferences between DMC and PIMC results, particularly 
for A^ > 26. PIMC chemical potentials show very promi- 
nent peaks at N=26, 29, 34 and 39, in contrast with the 
smooth behavior obtained with DMC. 

In this work we present new DMC calculations, im- 
proving our previous ones so as to get very precise 
results within our computational capacity. Specifically, 
we consider three aspects: the importance sampling func- 



tion, the time step adjustement and the estimation of the 
statistical errors. 

The DMC procedure is significantly improved when 
using a good importance sampling wave function, the 
main effect being the reduction of the variance of the 
stochastic procedure. We have used a Jastrow function 
with two- and three-body correlations: 
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Indices i,j,l run over the number of molecules in the 
cluster. This function is described in terms of five vari- 
ational parameterSj_25,pi, St, wt, and At- In our pre- 
vious calculations [10| we used the standard two-body 
Jastrow function <I>t = exp(u2). The present trial func- 
tion includes an enlarged the two-body variational space 
and also three-body correlations in the form suggested in 
Ref. [1^, which still has 0{N'^) computational complex- 
ity. 

DMC is based in a short-time approximation of 
the Green's function related to the imaginary time 
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Schrodinger equation. In this way, an initial wave func- 
tion ^xit — 0) evolves to the exact ground state wave 
function ^ at large t after many short-time steps r. We 
have used the 0{t^) approximation to the Green's func- 
tion as described in Refs. [13, [III i which provides ener- 
gies O(t^). The time step adjustment is the following: 
from calculations at the relative large steps 0.001 and 
0.0005K~^, we obtain the Richardson extrapolated value 



(4£;(0.0005) - £:(0.001)) 
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based on the r expansion E{t) = E{0) + Ct^ + ■ ■ ■ . This 
value turns out to be very close to the calculations with 
much smaller time steps, as it may be checked in the last 
three rows of Table |T1 This checks that the algorithm 
behaves as ©(r^), as expected, and suggests to use the 
value T = O.OOOIK"^ for massive calculations with a neg- 
ligible bias. 



TABLE I: Determination of optimal time step r from differ- 
ent evaluations of the binding energy B(N) of several clusters. 
The row labelled R.E. is the Richardson extrapolated value 
obtained from the previous two rows. The statistical standard 
deviation is indicated in parenthesis (error in the final digit 
shown). Energies and statistical errors arc in K. 

r B{10) B(20) 5(30) 

0.001 183.47(5) 559.28(17) 1006.4(3) 

0.0005 185.91(6) 566.56(17) 1020.0(4) 

R.E. 186.72(9) 568.99(28) 1024.5(5) 

0.0001 186.93(6) 569.16(12) 1025.2(2) 

0.00002 186.72(3) 569.48(7) 1024.8(1) 



A further improvement of the calculation regards the 
estimate of the statistical error. Because of the sequential 
Markov chain nature of Monte Carlo algorithms, succes- 
sive samples are strongly correlated, and the typical way 
of estimating the variance, cr^ = (-ff^) — (H)"^, may be 
too optimistic. To avoid these correlations we computed 
a number of times (10, typically), the binding energies, 
with independent and randomized runs, and estimate the 
variance from these results. This requires a considerable 
increasing in computational time, but the obtained stan- 
dard deviations are very precisely computed. Specifically, 
we have used 1000 walkers with 10^ steps plus 20000 sta- 
bilization steps in each walker. 

Hydrogen molecules interact through weak van der 
Waals forces that, nevertheless, are sufficiently strong to 
bound clusters with any number of molecules. Several 
forms have been derived to describe the p-H2"p-H2 inter- 
action. Two of them are of particular interest because 
they combine ab initio properties with properties of the 
gas (or solid) as well as experimental information from 
collisions, one due to Silvera and Goldman [l^, and the 
other to Buck et. al. [l^, hereafter referred to as SG 
and BHKOS, respectively. The main difference among 
them is that the former contains a repulsive long-range 
term (cg/r^) with the objective of providing an approxi- 



mation for the effective potential in a solid. Recent cal- 
culations have employed BHKOS potential SG po- 
tential or both [HI- We present here results 
with both interactions. 



TABLE II: DMC ground state (p-H2)jv binding energies (in 
K) obtained with BHKOS interaction [H]. 
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29 


978.50(10) 
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39 
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TABLE III: DMC ground state (p-H2)jv binding energies (in 
K) obtained with the SG interaction dSj. 
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The calculated DMC ground state binding energies are 
displayed in Tables [Hand [ml for BHKOS and SG poten- 
tials, respectively. As usual, the numbers in parenthesis 
are the errors in the final digit shown, and correspond 
to the standard deviation. The binding energies for the 
dimer have been obtained by numerical integration of the 
Schrodinger equation. 

The inclusion of triplet correlations in the importance 
sampling function leads to a noticeably improvement of 
the variational energies, as showed in Fig. [1] for BHKOS 
potential. The DMC energies are basically the same as 
those of Re f. llOl , with slightly more binding in the heav- 
ier clusters [20| . The main difference lies in the reduction 
of the standard deviation by typically a factor of 2. 

The total binding energies grow monotonically with 
the number of constituents. In order to determine an 
enhanced stability related to magic sizes it is convenient 
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FIG. 1: (Color online) Comparison of the binding energies 
per molecule obtained with diffusion Monte Carlo (DMC) and 
a variational Monte Carlo (VMC) calculation based on Eq. ((T)) 
with two- and three-body correlations. The interaction is 
BHKOS. The error bars are within the symbol size. 
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should be mentioned that even if the relative error of the 
total energies is around 10"'*, the relative error in the 
separation energies can be as high as 10^^, as a conse- 
quence of the strong cancelations appearing when com- 
puting II. So, even after our formidable numerical effort, 
the absolute error of /i for ~ 40 may be as high as 
0.5 K. Having this fact in mind, the only possible struc- 
ture, apart from N=13, is N=36 for BHKOS potential. 
However, /isg is ~ 1.70 ±0.6 K higher than its neighbors 
M35,37 ^'^d one cannot exclude that it could simply be 
a statistical fluctuation. Consequently the only magical 
cluster firmly established here is N=13, independent of 
the interaction. 



FIG. 3: (Color online) DMC and PIMC chemical potentials of 
(p-H2)jv clusters as a function of the number A'^ of molecules, 
calculated with SG interaction p_5j. Filled circles and squares 
are PIMC results of Ref. [ij, the solid line corresponds to 
PIMC results of Ref. [l3|]. Error bars have not been drawn 
for the clarity of presentation. 



to analyze the variation with N of the dissociation en- 
ergy or chemical potential, defined from the ground state 
energies E{N) as 



^jv = E{N ^1) - E{N) . 



(8) 



This quantity is plotted in Fig. [2] as a function of the 
number of molecules N. The main physical result of this 
figure is the presence of a neat peak at = 13, indicat- 
ing the magical character of this cluster. Although the 
two used interactions give different total energies, the 
BHKOS potential providing more binding than the SG 
one, this peak is present for both interactions. 
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FIG. 2: (Color online) DMC chemical potential (in K) of 
(p-H2)jv clusters as a function of the number N of molecules. 
The error bars are within the symbol size. 
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Beyond iV = 13 our calculations do not show any clear 
signal of local enhancement of the chemical potential. It 



In contrast, PIMC calculations show very prominent 
peaks at iV = 26, 29, 32, 34 and 39, as show in Fig. H 
The PIMC results at T = 0.5 K and 1.5 K are from 
Ref. [13, and those at T = 1 K from Ref. [H. Up 
to iV w 25, our DMC results are indistinguishable from 
the the PIMC ones of Ref. [li at T = 0.5 K or those of 
Ref. [p at T = 1 K. It is worth noticing that these PIMC 
results at r = 0.5 and 1 K are essentially identical in 
the calculated range of cluster sizes, with the noticeable 
exceptions of iV = 23, 35 and 36. However, it should be 
kept in mind that PIMC error bars have not been drawn 
in Fig. [3] for the sake of clarity. 

The existence of peaks in the chemical potential seems 
to be related to thermal effects. These could manifest in 
enhanced stability thresholds at finite temperature, sim- 
ilarly to what has been observed in ''He droplets [Uj. 
But according to Ref. [ll] such thermal effects should 
be associated to a coexistence of solid-like and liquid-like 
phases, with a dominance of the latter at low T, as a 
result of both the zero-point motion and quantum per- 
mutation exchanges. To this respect it is worth recalling 
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that while DMC may be affected by the constraint im- 
posed by the importance samphng function, PIMC has 
no such constraint. Actually, our importance sampling 
function is of the type used to describe liquid-like clus- 
ters. We have checked that the DMC ground state energy 
for the magical (p-H2)i3 cluster does not change employ- 
ing instead an importance sampling function where the 
molecules are localized at the vertex of the corresponding 
truncated MacKay icosahedra. This cluster is definitely 
liquid-like, in agreement with PIMC results. It would 



be interesting to perform a similar solid-like DMC calcu- 
lation for N > 26. Although very computationally de- 
manding, such a calculation could be useful to ascertain 
the phase of P-II2 clusters in this size region. 
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